num_metabolites = 6;
num_fluxes = 8;
kf = 10000.*rand(1, num_fluxes);
kr = 10000.*rand(1, num_fluxes);
epsilon = 0.0000001;
setpoint = 10;
Sf = [
        1 0 0 0 0 0 0 0;
        0 0 2 0 0 0 0 0;
        0 1 0 0 0 0 1 0;
        0 1 0 0 0 0 0 0;
        0 0 0 0 1 0 0 0;
        0 0 0 1 0 0 0 0
      ];
Sr = [
        0 1 1 0 0 0 0 0;
        0 0 0 1 1 0 0 0;
        0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 1 1;
        0 0 0 0 0 1 0 0;
        0 0 0 0 0 1 0 0
      ];
    
S = [
        1   -1	-1	0	0	0	0	0;
        0	0	2	-1	-1	0	0	0;
        0	1	0	0	-1	0	1	0;
        0	1	0	0	0	0	-1	-1;
        0	0	0	0	1	-1	0	0;
        0	0	0	1	0	-1	0	0
    ];
V = [3;2;1;1;1;1;-1;3];

[vf,vr,a]=steadystate(kf,kr,Sf,Sr,epsilon,setpoint,num_fluxes,num_metabolites)